Multivariate CAViaR

A More Insightful Approach to Risk Modeling

Abstract

This paper builds upon previous literature for modeling value-at-risk quantiles that uses an ARMA-style approach by adding ETFs as explanatory variables that are combined into principal component vectors at the forecast origin. By combining these principal component vectors along with various transformations of lagged autoregressive response variables, similar predictive accuracy is achieved (relative to the univariate model) and a clearer interpretation is gained. In fact, this thesis proposes an algorithm to detect changepoints in the explanatory variables that could be used as a means to measure changepoints in the economy. [DELETE IF THE CHANGEPOINT ARGUMENT DOESN’T MAKE SENSE.]

Background and Introduction

When modeling financial time series, simply considering the mean and the variance is insufficient. In fact, modeling a 1% or a 5% quantile of daily returns is a way to understand what happens on the worst days and to have a clearer picture of what might happen during a downturn. Indeed, finance theory suggests that a primary reason why the S&P500, which is a market-capitalization weighted index composed of the 500-largest publicly traded companies in the United States, has earned 6.8% inflation-adjusted pre-tax return with dividend reinvestment from January 1871 through April 2020 (CITE) is because of the risk of a significant downturn. Kerry Pechter at Forbes describes it as a premium for the fact that “stocks are riskier” and “more prone to price fluctuations in the short run” compared to lower risk investments. While the long-run picture of stock returns is a bit clearer, the short-run is a huge question - an investment manager using financial leverage to magnify returns (positive or negative) that could leave them in dire straits if their investments fell rapidly, despite a sound long-run strategy.

While there are other ways to understand and measure downside risk, a commonly accepted means is using value-at-risk (VaR). The metric is understood as follows: a one-day 1% VaR of -1 million dollars for a portfolio means that the portfolio will lose at least $1 million on the 1% worst trading days. A major advantage of VaR is that it distills a distribution of returns into one number. As such, VaR is often used in stress testing by regulatory agencies in the United States, the United Kingdom, and Europe (CITE).

Many of the approaches for modeling VaR rely on a semiparametric or a nonparametric historical simulation (ADD CITATION - BOUDOUKH). According to Robert Engle and Simone Manganelli in a 2004 paper, these methods are usually chosen for “empirical justifications rather than on sound statistical theory” (ADD CITATION - ENGLE). As such, they propose a framework called CAViaR that forecasts the VaR quantile directly using a conditional autoregressive quantile specification. This approach builds upon the statistical literature that extends linear quantile models to settings amenable to financial modeling, such as with heteroskedastic and nonstationary error distributions (ADD CITATION - KOENKER, PORTNOY).

Methods Used

A common problem in statistical analysis occurs when a training sample is significantly different than that of the test sample. Initial motivations for this paper involved analyzing two stocks - Amazon (ticker: AMZN) and Proctor & Gamble (ticker: PG) and their performance during the great recession. A relevant question of a financial institution would understandably be how their risk model performed during 2008, a highly volatile period which was driven by the “worst financial crisis since the Great Depression” according to Gary Becker (CITE), a Nobel-prize winning Economist. Interestingly, the forecast for Amazon was fairly accurate whereas the forecast for PG was not. One reason for this could be the fact that a stock like Amazon was highly volatile during the training sample, which included data from previous years (HOW MANY?), but PG was fairly calm. How would it be possible for a univariate model such as CAViaR, that does not explictly account for other factors, to forecast well? What if a volatile stock such as AMZN was included into the forecast for PG - would it improve the forecast?

Thus, the idea of combining stocks into a multivariate setting to capture correlations and better forecast risk was formed. A natural choice appeared to be the diffusion index model, originally developed by Stock and Watson for predicting conditional means [TWO CITATIONS]. The model for forecasting the conditional mean is specified below.

Diffusion Index Model

A useful means of predicting stock movements in the future is the Stock and Watson diffusion index. The model is outlined below, which is adapted from Multivariate Time Series Analysis With R and Financial Applications by Ruey S. Tsay [ADD CITATION].

There are two relevant equations:

\[ \boldsymbol{z_t} = \boldsymbol{Lf_t} + \boldsymbol{\epsilon_t} \]

and

\[ y_{t+h} = \boldsymbol{\beta^\prime f_t} + \boldsymbol{e_{t+h}} \]

In the first equation \(\boldsymbol{z_t} = (z_{1t}, ...., z_{kt})^\prime\) is an observed time series with mean 0, \(\boldsymbol{f_t}\) is an m-dimensional vector of common factors with mean 0 and identity covariance matrix, \(\boldsymbol{L}\) is a \(k\) x \(m\) loading matrix, and \(\boldsymbol{\epsilon_t}\) is a i.i.d. sequence of random vectors with mean 0 and covariance matrix \(\boldsymbol{\Sigma_e}\).

In the second equation, which represents the h-step ahead prediction based on \(\boldsymbol{f_t}\), \(y_t\) is the scalar time series of interest, \(h\) is the forecast horizon, \(\boldsymbol{\beta}\) represents the vector of coefficients, and \(e_t\) is a sequence of uncorrelated random variables with mean 0 and constant variance.

To model the data, principal component analysis is performed on the covariates described below to obtain an estimate of \(\boldsymbol{f_t}\). Then the \(\boldsymbol{\beta}\) coefficients are estimated using ordinary least squares.

CAViaR Model

However, work needed to be done to align the diffusion index model with the CAViaR model, which is defined below. The following variables are required for use in the CAViaR model:

A vector of portfolio returns

Adaptive CAViaR Model

\[ f_t(\beta_1) = f_{t-1}(\beta_1) + \beta_1\left[\left(1+ \exp(G[y_{t-1} - f_{t-1}(\beta_1)]) \right)^{-1} - \theta \right] \]

Symmetric Absolute Value CAViaR Model

\[ f_t(\boldsymbol{\beta}) = \beta_1 + \beta_2f_{t-1}(\boldsymbol{\beta}) + \beta_3|y_{t-1}| \]

Asymmetric Slope CAViaR Model

\[ f_t(\boldsymbol{\beta}) = \beta_1 + \beta_2f_{t-1}(\boldsymbol{\beta}) + \beta_3(y_{t-1})^+ + \beta_4(y_{t-1})^- \]

Indirect GARCH (1,1) CAViaR Model

\[ f_t(\boldsymbol{\beta}) = (\beta_1 + \beta_2f_{t-1}^2(\boldsymbol{\beta}) + \beta_3y_{t-1}^2)^{1/2} \]

[ADD SPECIFICATION]

As such, the following model is used to combine the two ideas:

VAR = INTERCEPT + PAST VALUES + AR COMPONENT + ERROR

Data Used

The response variable used in this analysis is SPY, which is an exchange-traded fund that aims to track the performance of the S&P 500, which is discussed above. It is broadly used as a bellweather of the U.S. economy (CITE), and has the advantage of avoiding survivorship bias - while an individual stock might go bankrupt or merge with another, it’s reasonable to assume that these issues do not apply with an ETF.

Following this logic, there are several classes of response variables used in this analysis. The first group is a set of U.S. sector ETFs obtained from Seeking Alpha using the link below. As with the response variable, these ETFs were publicly traded throughout the Great Recesssion.

https://seekingalpha.com/etfs-and-funds/etf-tables/sectors

1. Utilities (XLU)
2. Consumer Staples (XLP)
3. Healthcare (XLV)
4. Technology (XLK)
5. Consumer Discretionary (XLY)
6. Industrial (XLI)
7. Financial Services (XLF)
8. Basic Materials (XLB)
9. Energy (XLE)

The second is Global Sector ETFs, also obtained from Seeking Alpha. The rationale for including these is that perhaps some global exposure is useful in understanding the broader market.

1. Utilities (JXI)
2. Consumer Staples (KXI)
3. Healthcare (IXJ)
4. Telecommunications (IXP)
5. Technology (IXN)
6. Consumer Discretionary (RXI)
7. Industrial (EXI)
8. Financial Services (IXG)
9. Basic Materials (MXI)
10. Energy (IXC)

The third is bond ETFs, also obtained from [CITE]. Like the previous two, these ETFs potentially contain forward-looking information about the stock market. Link below:

1. iShares 1-3 Year Treasury Bond Fund (SHY)
2. iShares 7-10 Year Treasury Bond Fund (IEF)
3. iShares 20+ Year Treasury Bond Fund (TLT)
4. iShares iBoxx $ Investment Grade Corporate Bond ETF (LQD)

Lastly, all of the above are run together.

Results

The working directory was changed to /Users/stevenmoen/Documents/GitHub/CAViaR inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

Big Simulation Function

Problem Solving on 4.29.2020

2008 Ending

U.S. ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.046	1.1	1.371	1.265	0.208	0.213	0.219	0.355
VaR Breaks (%)	0.200	0.2	0.212	0.208	0.028	0.028	0.028	0.060
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.319	1.344	1.768	1.385	0.651	0.654	0.640	0.956
VaR Breaks (%)	0.260	0.236	0.340	0.260	0.076	0.076	0.064	0.160
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.623	1.564	1.738	1.534	1.077	1.066	1.068	1.366
VaR Breaks (%)	0.328	0.312	0.348	0.332	0.144	0.156	0.140	0.224
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Global ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.740	0.761	0.841	0.867	0.208	0.213	0.219	0.355
VaR Breaks (%)	0.108	0.112	0.120	0.108	0.028	0.028	0.028	0.060
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.160	1.173	2.157	1.283	0.651	0.654	0.640	0.956
VaR Breaks (%)	0.176	0.172	0.412	0.184	0.076	0.076	0.064	0.160
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.517	1.517	1.791	1.523	1.077	1.066	1.068	1.366
VaR Breaks (%)	0.284	0.284	0.348	0.288	0.144	0.156	0.140	0.224
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Commodity ETFs

Bond ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.891	1.017	1.271	1.090	0.208	0.213	0.219	0.355
VaR Breaks (%)	0.156	0.168	0.224	0.128	0.028	0.028	0.028	0.060
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.338	1.316	1.659	1.325	0.651	0.654	0.640	0.956
VaR Breaks (%)	0.236	0.236	0.316	0.256	0.076	0.076	0.064	0.160
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.641	1.575	1.840	1.586	1.077	1.066	1.068	1.366
VaR Breaks (%)	0.308	0.304	0.364	0.308	0.144	0.156	0.140	0.224
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

All ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.736	0.737	1.733	0.863	0.208	0.213	0.219	0.355
VaR Breaks (%)	0.104	0.108	0.216	0.104	0.028	0.028	0.028	0.060
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.148	1.236	2.391	1.371	0.651	0.654	0.640	0.956
VaR Breaks (%)	0.168	0.208	0.400	0.220	0.076	0.076	0.064	0.160
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	1.521	1.549	1.797	1.644	1.077	1.066	1.068	1.366
VaR Breaks (%)	0.284	0.288	0.344	0.292	0.144	0.156	0.140	0.224
Note:
Calculated using 250 trading days from 2008-01-04 to 2008-12-30

2010 Ending

U.S. ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.128	0.119	0.155	0.128	0.079	0.08	0.086	0.191
VaR Breaks (%)	0.000	0.000	0.044	0.000	0.020	0.02	0.016	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.369	0.369	0.468	0.379	0.336	0.336	0.343	0.492
VaR Breaks (%)	0.028	0.028	0.084	0.028	0.052	0.052	0.048	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.56	0.561	0.710	0.572	0.547	0.549	0.546	0.690
VaR Breaks (%)	0.08	0.072	0.132	0.080	0.080	0.088	0.084	0.028
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Global ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.133	0.136	0.110	0.127	0.079	0.08	0.086	0.191
VaR Breaks (%)	0.000	0.000	0.004	0.000	0.020	0.02	0.016	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.386	0.371	0.449	0.401	0.336	0.336	0.343	0.492
VaR Breaks (%)	0.024	0.024	0.068	0.024	0.052	0.052	0.048	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.606	0.606	0.664	0.590	0.547	0.549	0.546	0.690
VaR Breaks (%)	0.040	0.040	0.128	0.068	0.080	0.088	0.084	0.028
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Commodity ETFs

Bond ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.131	0.13	0.129	0.128	0.079	0.08	0.086	0.191
VaR Breaks (%)	0.000	0.00	0.016	0.004	0.020	0.02	0.016	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.373	0.372	0.494	0.378	0.336	0.336	0.343	0.492
VaR Breaks (%)	0.028	0.028	0.104	0.028	0.052	0.052	0.048	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.566	0.565	0.659	0.578	0.547	0.549	0.546	0.690
VaR Breaks (%)	0.076	0.076	0.128	0.080	0.080	0.088	0.084	0.028
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

All ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.146	0.136	0.114	0.128	0.079	0.08	0.086	0.191
VaR Breaks (%)	0.000	0.000	0.000	0.000	0.020	0.02	0.016	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.394	0.394	0.436	0.397	0.336	0.336	0.343	0.492
VaR Breaks (%)	0.024	0.024	0.068	0.028	0.052	0.052	0.048	0.000
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.595	0.594	0.664	0.596	0.547	0.549	0.546	0.690
VaR Breaks (%)	0.044	0.044	0.120	0.076	0.080	0.088	0.084	0.028
Note:
Calculated using 250 trading days from 2010-01-05 to 2010-12-30

2014 Ending

U.S. ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.083	0.075	0.173	0.071	0.061	0.057	0.063	0.061
VaR Breaks (%)	0.000	0.000	0.072	0.004	0.008	0.004	0.012	0.028
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.242	0.241	0.229	0.238	0.226	0.218	0.225	0.240
VaR Breaks (%)	0.024	0.028	0.044	0.032	0.052	0.048	0.052	0.056
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.388	0.362	0.368	0.367	0.367	0.359	0.364	0.368
VaR Breaks (%)	0.060	0.056	0.080	0.076	0.116	0.104	0.112	0.132
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Global ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.076	0.074	0.158	0.071	0.061	0.057	0.063	0.061
VaR Breaks (%)	0.000	0.000	0.096	0.000	0.008	0.004	0.012	0.028
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.240	0.245	0.241	0.245	0.226	0.218	0.225	0.240
VaR Breaks (%)	0.024	0.028	0.040	0.024	0.052	0.048	0.052	0.056
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.377	0.358	0.388	0.358	0.367	0.359	0.364	0.368
VaR Breaks (%)	0.056	0.048	0.084	0.068	0.116	0.104	0.112	0.132
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Commodity ETFs

Bond ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.079	0.075	0.072	0.070	0.061	0.057	0.063	0.061
VaR Breaks (%)	0.000	0.000	0.000	0.004	0.008	0.004	0.012	0.028
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.241	0.237	0.231	0.237	0.226	0.218	0.225	0.240
VaR Breaks (%)	0.028	0.012	0.040	0.024	0.052	0.048	0.052	0.056
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.370	0.364	0.371	0.352	0.367	0.359	0.364	0.368
VaR Breaks (%)	0.056	0.044	0.072	0.064	0.116	0.104	0.112	0.132
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

All ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.087	0.073	0.107	0.079	0.061	0.057	0.063	0.061
VaR Breaks (%)	0.008	0.008	0.052	0.008	0.008	0.004	0.012	0.028
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.241	0.256	0.320	0.246	0.226	0.218	0.225	0.240
VaR Breaks (%)	0.024	0.032	0.084	0.028	0.052	0.048	0.052	0.056
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.371	0.359	0.370	0.361	0.367	0.359	0.364	0.368
VaR Breaks (%)	0.056	0.044	0.072	0.056	0.116	0.104	0.112	0.132
Note:
Calculated using 250 trading days from 2014-01-03 to 2014-12-30

2016 Ending

U.S. ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.087	0.093	0.085	0.09	0.078	0.082	0.078	0.077
VaR Breaks (%)	0.020	0.020	0.028	0.02	0.012	0.020	0.012	0.004
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.278	0.297	0.257	0.290	0.238	0.238	0.234	0.264
VaR Breaks (%)	0.064	0.056	0.080	0.068	0.032	0.040	0.028	0.032
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.419	0.415	0.401	0.429	0.370	0.373	0.368	0.414
VaR Breaks (%)	0.104	0.100	0.120	0.108	0.088	0.092	0.096	0.072
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Global ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.09	0.098	0.085	0.087	0.078	0.082	0.078	0.077
VaR Breaks (%)	0.02	0.020	0.028	0.016	0.012	0.020	0.012	0.004
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.279	0.261	0.253	0.291	0.238	0.238	0.234	0.264
VaR Breaks (%)	0.068	0.040	0.048	0.080	0.032	0.040	0.028	0.032
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.419	0.400	0.400	0.442	0.370	0.373	0.368	0.414
VaR Breaks (%)	0.104	0.096	0.116	0.112	0.088	0.092	0.096	0.072
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Commodity ETFs

Bond ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.076	0.089	0.108	0.085	0.078	0.082	0.078	0.077
VaR Breaks (%)	0.012	0.024	0.028	0.016	0.012	0.020	0.012	0.004
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.257	0.265	0.273	0.273	0.238	0.238	0.234	0.264
VaR Breaks (%)	0.040	0.044	0.060	0.060	0.032	0.040	0.028	0.032
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.402	0.395	0.397	0.415	0.370	0.373	0.368	0.414
VaR Breaks (%)	0.100	0.096	0.108	0.084	0.088	0.092	0.096	0.072
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

All ETFs

Comparison of VaR Methods for a 1% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.091	0.097	0.092	0.083	0.078	0.082	0.078	0.077
VaR Breaks (%)	0.020	0.020	0.028	0.016	0.012	0.020	0.012	0.004
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 5% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.26	0.265	0.283	0.272	0.238	0.238	0.234	0.264
VaR Breaks (%)	0.06	0.048	0.092	0.056	0.032	0.040	0.028	0.032
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Comparison of VaR Methods for a 10% VaR

	M1	M2	M3	M4	M5	M6	M7	M8
Losses	0.425	0.397	0.401	0.441	0.370	0.373	0.368	0.414
VaR Breaks (%)	0.096	0.096	0.112	0.100	0.088	0.092	0.096	0.072
Note:
Calculated using 250 trading days from 2016-01-06 to 2016-12-30

Conclusions and Future Work

The problem of how to predict a low quantile of a stock’s log return when the training sample is substantially different than the test scenario is an enormously difficult problem. Almost axiomatically, the distribution is nonstationary over time. How is it possible to predict the return of an index like the S&P500 during a period of market turmoil such as the great recession? While the CAViaR model performs comparatively well during times of stress, it performs about the same as the multivariate model during more benign economic periods.

This conclusion drawn from the above results might support the notion of combining the two models in some sort of a mixture model - aiming to use the basket of ETFs during good times, and use the CAViaR ARMA specification during bad times. The approach of using ETFs allows a prediction based on forward-looking expectations of fundamental factors. Indeed, ETFs are just baskets of individual stocks or bonds, and those securities are (in theory) based on rational expectations about future resources, market conditions, etc - the microfoundations of what drives our economy. The ARMA specification, while practically and statistically sound, is contradicted by economic theory and practice - the weak form of the efficient market hypothesis states that it is impossible to forecast future values of asset prices using past values. But perhaps this view is incomplete.

Any model that attempts to capture relationships in the real world will only work until an omitted variable is found. The elegance of the multivariate CAViaR model is that it provides insight into why a prediction is wrong; the change in the angle between resultant vectors is a sensible measurement of economic changepoints (CHANGE). However, errors in the world are costly, and it’s wishful thinking to say that explaining why the error occured is sufficient.

As such, for future work it’s worth exploring the notion of weighting an ARMA-approach more heavily when predictions using fundamentals were too high, then not only would this after-the-fact recognition be achieves, but also a hierarchical model that captures fundamental relationships in the economy and potentially changes our understanding of asset prices in general - a synthesis between Keynes’ animal spirits during a time of severe crisis; where a model cannot explain shifts, and a more rational world that explains other periods. In addition to significant predictive power because of the switching between the two worlds, there’s also an elegant explanation; a way to explain changes in the usefulness of the underpinnings in the economy. Because of the flexibility of the model, it’s entirely possible that a whole gamut of variables could be tossed in and backtested to when “changepoints” occurred.

Additional future work involves developing theoretical guarantees on the parameters in the multivariate CAViaR model. One advantage of both the diffusion index model and the CAViaR model is that both have theorems about asymptotic normality and consistency [CONFIRM].

Code Appendix

knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/',
                      echo=FALSE, warning=FALSE, message=FALSE, cache = TRUE)
# Read in relevant libraries
library(microbenchmark)
library(data.table)
library(quantmod)
library(ggplot2)
library(tseries)
library(zoo)
library(magrittr)
library(dplyr)
library(kableExtra)
library(formattable)
library(quantreg)
library(MTS)
library(plot3D)
library(citr)

# Set up working directory
# setwd("~/Documents/GitHub/CaviaR")

# source('caviar_SM.R')
source('~/Documents/GitHub/CaviaR/caviar_SM.R')
# This code below is for use in the CAViaR sections.
# Here is code that I'll wrap some parts in to avoid superfluous output
quiet <- function(x) { 
  sink(tempfile()) 
  on.exit(sink()) 
  invisible(force(x)) 
} 
#' This is a function which pulls data for use in the CAViaR model
#'
#' @param symbol - symbol to pull
#' @param compl_case - defaults to true...only includes complete cases in the data
#' @param adj_close - use adjusted closing prices. Default is yes.
#' @param log_return - use log return? Default is yes.
#'
#' @return - a data frame which can be fed into later functions
#' @export
#'
#' @examples - data_pull("SPY")
data_pull = function(symbol, compl_case = 1, adj_close = 1, log_return = 1, start_date = "1900-01-01", end_date = Sys.Date()){
  # Pull in data from quantmod
  response_pull = getSymbols(symbol, auto.assign = FALSE, from = start_date, to = end_date)
  # Get adjusted closing price
  if (adj_close == TRUE){
    df = Ad(response_pull)
  } else {
    df = Cl(response_pull)
  }
  # Return complete cases only 
  if (compl_case == TRUE){
    df = df[complete.cases(df), ]
  } else{
    df = df
  }
  # Calculate log return of data
  if (log_return == TRUE){
    lr = log(df[,1]/shift(df[,1], 1, type = "lag"))
    # Combine data
    df_out = cbind(df, lr)
    # Rename the data 
    colnames(df_out) <- c(sym=symbol, paste0(symbol, "_log_return"))
  } else{
    df_out = df
  }
  # Return data
  return(df_out)
}


#' Pull the data and run the CAViaR function on it
#'
#' @param input_data - data to use in the function
#' @param range_data - range of the data to use
#'
#' @return - a list of values from the caviar function
#' @export
#'
#' @examples - caviar_pull(spy)
caviar_pull = function(input_data, range_data = (2:dim(input_data)[1])){
  # Run the caviar data
  caviar <- caviarOptim(input_data[range_data,2])
  return(caviar)
}

#' Function for producing rolling predictions
#' Model 1 = Symmetric Absolute Value, 2 = Asymmetric slope, 3 = Indirect GARCH, 4 = Adaptive
#'
#' @param input_data - input data from the previous function
#' @param range_data - range of the data to consider
#' @param nfcst - number of forecasts to make
#' @param model - model to use (integers 1 through 4). Defaults to 1. 
#' @param level - level of significance to use.
#' @param G - argument for the k parameter in the 4th model (adaptive). Default is 5
#'
#' @return - an xts object which contains rolling CAViaR predictions
#' @export
#'
#' @examples - rolling_predictions(spy, nfcst = 22)
rolling_predictions = function(input_data, range_data = (2:dim(input_data)[1]), nfcst = 250, model =1, level = 0.01, G = 5, col = 2){
  # Run the varpredict function
  varpredict <- rollapplyr(input_data[range_data,col], length(range_data) - nfcst, caviarOptim, model, level, predict = 1, k = G) %>% lag
  # Eliminate NAs
  # pred_no_na = na.omit(varpredict)
  # Return the data
  # return(pred_no_na)
  return(varpredict)
}

#' Function to Calculate Loss from the above predictions
#'
#' @param symbol - symbol to work with from quantmod. Must be in quotations to work
#' @param start_dt - start date of the data to build the forecast on 
#' @param end_dt - end date of the data to build the forecast on  
#' @param nfcst - number of data points to use in the forecast
#' @param model - model to use. Defaults to 1
#' @param level - level of significance. Defaults to 1%
#' @param G - argument for the k parameter in the 4th model (adaptive). Default is 5
#'
#' @return - loss using absolute value
#' @export - a plot of the data
#'
#' @examples
loss_calc_uv = function(symbol, start_dt, end_dt, nfcst, model = 1, level = 0.01, G = 5){
  # Pull in the data
  raw_data = data_pull(symbol, start_date = start_dt, end_date = end_dt)
  # Forecast based on the data
  fcst = na.omit(rolling_predictions(raw_data, nfcst = nfcst, model = model, level = level, G = G))*(-1)
  # Extract actuals
  act = tail(raw_data, n = nfcst)[,2]
  # Join the two together and rename
  join = merge(fcst,act,all=TRUE)
  colnames(join) <- c("Fcst_VaR", "Act_Return")
  # print(join)
  # Calculate the losses
  loss = abs(sum(ifelse(act > fcst, level, (-1)*(1-level))))
  # Plot the data
  plot = plot.xts(join, col = c("red", "black"), lty = c(2,1), main = "Log Return from the SPY vs. Fcst. VaR",grid.col = NA, legend.loc = "bottomleft")
  return(list(loss, plot, act, fcst))
}

#' This is a function which creates a data frame for the response and explanatory variables that we'll feed into the diffusion index
#'
#' @param symbol_list - a list of symbols recognizable by the 
#' @param resp_var - the response variable we'd like to forecast; default is SPY
#' @param compl_case - defaults to true...only includes complete cases in the data
#' @param adj_close - use adjusted closing prices for the explanatory variables? default is 1 for YES
#' @param resp_adj_close - use adjusted closing prices for the explanatory variables? default is 1 for YES
#' @param start_date - starting data to use
#' @param end_date - ending date of the data
#' @param lag_pred - do we lag the predictions? It is STRONGLY recommended that this is 0
#'
#' @return - a data frame which can be fed into the SWfore function
#' @export
#'
#' @examples - diff_index_df(c("XLF", "XLE", "PSCT", "XLV", "VPU", "XLP", "IGF", "XWEB", "PPTY"))
diff_index_df = function(symbol_list, resp_var = "SPY", compl_case = 1, adj_close = 1, resp_adj_close = 1, start_date = "1900-01-01", end_date = Sys.Date(), lag_pred = 1){
  # Pull in response variable
  response_pull = getSymbols(resp_var, auto.assign = FALSE, from = start_date, to = end_date)
  # Get adjusted closing price
  if (resp_adj_close == TRUE){
    diff_df = Ad(response_pull)
  } else {
    diff_df = Cl(response_pull)
  }
  # Loop through the symbols and join in data
  for (i in 1:length(symbol_list)){
    # Pull closing price
    expl_pull = getSymbols(symbol_list[i], auto.assign = FALSE, from = start_date, to = end_date)
    # Extract closing price - 4th element
    if (adj_close == TRUE){
      expl_cl = Ad(expl_pull)
    } else {
      expl_cl = Cl(expl_pull)
    }
    # New code for 4.16.2020 - lag the explanatory variables
    if (lag_pred == TRUE){
      # Lag the explanatory variables by 1
      lag_exp = lag(expl_cl, 1)
      # Append the first lag to the data frame
      diff_df = merge(diff_df, lag_exp, join = "left", fill = NA)
    } else{
      # Return the data frame without lags
      diff_df = merge(diff_df, expl_cl, join = "left", fill = NA)
    }
  }
  if (lag_pred == TRUE){
    # Chop off the first row
    diff_df = diff_df[-1,]
  }
  else {
    print("PLEASE NOTE - the explanatory variables in this DF are NOT lagged. Be careful to avoid look-ahead bias!")
  }
  # Return complete cases only 
  if (compl_case == TRUE){
    diff_df_out = diff_df[complete.cases(diff_df), ]
  } else{
    diff_df_out = diff_df
  }
  
  return(diff_df_out)
}

#' Converts a diff_df into a data frame with approximate percentage changes diff(log(diff_df))
#'
#' @param diff_df - output of the diff_index_df function with complete cases
#'
#' @return - retuns the differenced data
#' @export
#'
#' @examples - pc_diff_index(test_compl) 

pc_diff_index = function(diff_df){
  # Difference the log of the data
  pc_diff_index = diff(log(diff_df))
  # Remove the first row
  pc_diff_index_out = pc_diff_index[-1,]
  return(pc_diff_index_out)
}

#' Below is the modified diffusion index code.
#'
#' @param y - response variable
#' @param x - predictor variables
#' @param orig - forecast origin
#' @param m - number of diffusion indexes used
#' @param tau - VaR level to use; must be between 0 and 1
#' @param end - specifies an alternate ending value
#' @param print_mdl - print the model summary and the MSE
#'
#' @return - returns a list of variables for use in the diffusion index
#' @export
#'
#' @examples
mod_di = function (y, x, orig, m, tau, end = NULL, print_mdl = 0) 
{
  # Converts the response variables into a matrix
  if (!is.matrix(x)) 
      x = as.matrix(x)
  # nT is number of t time-steps
  nT = dim(x)[1]
  # Add a line to establish the number of data points used in the test.
  if (is.null(end) != TRUE){
    nT = end
  }
  # k is the number of diffusion indices used
  k = dim(x)[2]
  # Sanity checks to ensure that the origin isn't past the number of time points
  if (orig > nT) 
      orig = nT
  # Makes sure that there aren't more predictors than there variables in the dataset
  if (m > k) 
      m = k
  # Makes sure there are at least some variables
  if (m < 1) 
      m = 1
  # Subdivides the dataframe
  x1 = x[1:orig, ]
  # Calculates means of each row
  me = apply(x1, 2, mean)
  # Calculates standard deviations of each column
  se = sqrt(apply(x1, 2, var))
  # Creates a matrix x1, which normalizes all the columns. 
  # This may be an issue since it assumes that the distribution is sufficiently described by the first two moments
  x1 = x
  for (i in 1:k) {
      x1[, i] = (x1[, i] - me[i])/se[i]
  }
  V1 = cov(x1[1:orig, ])
  # Performs an eigen decomposition
  m1 = eigen(V1)
  # Selects eigenvalues
  sdev = m1$values
  # Selects eigenvectors
  M = m1$vectors
  # Makes a smaller matrix
  M1 = M[, 1:m]
  # This is the diffusion index model - [orig x p]*[p x m] = [orig x m]
  Dindex = x1 %*% M1
  # Cut down both the response and predictors to be a reasonable size
  y1 = y[1:orig]
  DF = Dindex[1:orig, ]
  # Apply the linear model - HERE is the key.
  # mm = lm(y1 ~ DF) - old function
  mm = rq(y1 ~ DF, tau = tau)
  # Print the data
  if (print_mdl == 1){
    print(summary(mm))
  }
  # Puts coefficients in a matrix
  coef = matrix(mm$coefficients, (m + 1), 1)
  # Initializes yhat variables and MSE
  yhat = NULL
  MSE = NULL
  if (orig < nT) {
    # Creates a nfcst by (m+1) matrix
    newx = cbind(rep(1, (nT - orig)), Dindex[(orig + 1):nT, 
        ])
    # [nfcstx(m+1)]*[(m+1)x1] = [nfcstx1]
    yhat = newx %*% coef
    # Calculates errors
    err = y[(orig + 1):nT] - yhat
    MSE = mean(err^2)
    if (print_mdl == 1){
      cat("MSE of out-of-sample forecasts: ", MSE, "\n")
    }
  }
  SWfore <- list(coef = coef, yhat = yhat, MSE = MSE, loadings = M1, 
      DFindex = Dindex)
}

#' Below is the modified diffusion index code to include lagged variables.
#'
#' @param y - response variable
#' @param x - predictor variables
#' @param orig - forecast origin
#' @param m - number of diffusion indexes used
#' @param tau - VaR level to use; must be between 0 and 1
#' @param ar_tf - AR transformation type. (1 - no transformation,
#' 2 - absolute value, 3 - asymmetric slope)
#' @param p - number of AR lags to include. Default is one.
#' @param print_mdl - option to print the model summary to make sure everytning is ok. 0 is default.
#' @param model - model type (1 - SAV, 2 - AS, 3 - GARCH, 4 - ADAPTIVE) 
#'
#' @return - returns a list of variables for use in the diffusion index
#' @export
#'
#' @examples
mod_di_wl = function (y, x, orig, m, tau, ar_tf = 1, p = 1, print_mdl = 0, model = 1, end = NULL) 
{
  # Converts the response variables into a matrix
  if (!is.matrix(x)) 
      x = as.matrix(x)
  # nT is number of t time-steps
  nT = dim(x)[1]
  # Add a line to establish the number of data points used in the test.
  if (is.null(end) != TRUE){
    nT = end
  }
  # k is the number of diffusion indices used
  k = dim(x)[2]
  # Sanity checks to ensure that the origin isn't past the number of time points
  if (orig > nT) 
      orig = nT
  # Makes sure that there aren't more predictors than there variables in the dataset
  if (m > k) 
      m = k
  # Makes sure there are at least some variables
  if (m < 1) 
      m = 1
  # Subdivides the dataframe
  x1 = x[1:orig, ]
  # Calculates means of each row
  me = apply(x1, 2, mean)
  # Calculates standard deviations of each column
  se = sqrt(apply(x1, 2, var))
  # Creates a matrix x1, which normalizes all the columns. 
  # This may be an issue since it assumes that the distribution is sufficiently described by the first two moments
  x1 = x
  for (i in 1:k) {
      x1[, i] = (x1[, i] - me[i])/se[i]
  }
  V1 = cov(x1[1:orig, ])
  # Performs an eigen decomposition
  m1 = eigen(V1)
  # Selects eigenvalues
  sdev = m1$values
  # Selects eigenvectors
  M = m1$vectors
  # Makes a smaller matrix
  M1 = M[, 1:m]
  # This is the diffusion index model - [orig x p]*[p x m] = [orig x m]
  Dindex = x1 %*% M1
  # Cut down both the response and predictors to be a reasonable size
  y1 = y[1:orig]
  DF = Dindex[1:orig, ]
  # Copy the data frame
  DF_wl = Dindex
  # Lag the y-variable
  for (i in 1:p){
    # Create a lagged variable
    lag_var = lag(y, i)
    # Append the first lag to the data frame
    DF_wl = cbind(DF_wl,lag_var)
  }
  # Identify the right columns
  l_ar = ncol(DF_wl)
  f_ar = l_ar - p + 1
  # Keep the last columns kept to the side
  all_lag = DF_wl[,(f_ar:l_ar)]
  # Cut off the first row to avoid NA's
  DF_trim = DF_wl[1:orig,]
  # Rename the columns
  # Here's the new function with an untransformed AR(p) lag
  if (ar_tf == 1){
    # Incorporate everything in to an input data frame
    df_in = cbind(y1[-(1:p)], DF_trim[-(1:p),])
    # Rename the columns
    # Initialize a character vector
    nvec = c(rep(0, 1+m+p))
    # Populate the vector - first value is the response
    nvec[1] <- names(y)
    # Next are the diffusion indices
    for (i in 1:m){
      nvec[i+1] = paste0("Diff_Index_", i)
    }
    # Next are the lagged variables
    for (i in 1:p){
      nvec[i+1+m] = paste0("Lag_", i)
    }
    # Assign the names
    names(df_in) <- nvec
    # Run the model
    mm = rq(df_in[,1] ~ df_in[,-1], tau = tau)
  }
  # Here's the new function with an SAV AR(p) lag
  if (ar_tf == 2){
    # Incorporate everything in to an input data frame
    df_in = cbind(y1[-(1:p)], DF_trim[-(1:p),-(f_ar:l_ar)], abs(DF_trim[-(1:p),(f_ar:l_ar)]))
    # Rename the columns
    # Initialize a character vector
    nvec = c(rep(0, 1+m+p))
    # Populate the vector - first value is the response
    nvec[1] <- names(y)
    # Next are the diffusion indices
    for (i in 1:m){
      nvec[i+1] = paste0("Diff_Index_", i)
    }
    # Next are the lagged variables
    for (i in 1:p){
      nvec[i+1+m] = paste0("Lag_", i)
    }
    # Assign the names. Note that this is a matrix
    names(df_in) <- nvec
    # Run the model
    mm = rq(df_in[,1] ~ df_in[,-1], tau = tau)
  }
  # Here's the new function with an asymmetric slope for the AR(1) lag
  # Indicator; 0 if percent change is negative, 1 if it's positive
  # indi = ifelse(DF_trim[,ar] < 0, 0, 1)
  if (ar_tf == 3){
    # Create a matrix of indicators
    indi_mat = matrix(0, nrow(DF_wl), p)
    # Generalize the above code
    for (i in 1:p){
      # Populate the indicator
      indi_mat[,i] = ifelse(DF_wl[,f_ar + i - 1] < 0, 0, 1)
    }
  }
  # Fitting the regression
  if (ar_tf == 3){
    # Incorporate everything in to an input data frame
    df_in = cbind(y1[-(1:p)], DF_trim[-(1:p),-(f_ar:l_ar)], DF_trim[-(1:p),(f_ar:l_ar)], indi_mat[((p+1):orig),])
    # Rename the columns
    # Initialize a character vector
    nvec = c(rep(0, 1+m+2*p))
    # Populate the vector - first value is the response
    nvec[1] <- names(y)
    # Next are the diffusion indices
    for (i in 1:m){
      nvec[i+1] = paste0("Diff_Index_", i)
    }
    # Next are the lagged variables
    for (i in 1:p){
      nvec[i+1+m] = paste0("Lag_", i)
    }
    # Last are the positive indicator variables
    for (i in 1:p){
      nvec[i+1+m+p] = paste0("Pos_Val_for_Lag_", i)
    }
    # Assign the names. Note that this is a matrix
    names(df_in) <- nvec
    # Run the model
    mm = rq(df_in[,1] ~ df_in[,-1], tau = tau)
    # mm = rq(y1[-(1:p)] ~ DF_trim[-(1:p),-(f_ar:l_ar)] + DF_trim[-(1:p),(f_ar:l_ar)] + indi_mat[((p+1):orig),], tau = tau)
    # Add a different line to account for the indicator variable
    # intercept + m + 2*nlag to account for the number of indicator variables
    coef = matrix(mm$coefficients, (1 + m + 2*p), 1)
  }
  if (print_mdl == 1){
    print(summary(mm))
  }
  # Puts coefficients in a matrix - added the AR terms
  # coef = matrix(mm$coefficients, (m + 1), 1)
  if (ar_tf != 3){
    coef = matrix(mm$coefficients, (1 + m + p), 1)
  }
  # Initializes yhat variables and MSE
  yhat = NULL
  loss = NULL
  if (orig < nT) {
    # Creates a nfcst by (m+2) matrix
    # Add on the lagged variables
    newx = cbind(rep(1, (nT - orig)), Dindex[(orig + 1):nT, ], all_lag[(orig+1):nT,])
    # Incorporate lagged variables
    if (ar_tf == 3){
      newx = cbind(rep(1, (nT - orig)), Dindex[(orig + 1):nT, ], all_lag[(orig+1):nT,], indi_mat[(orig+1):nT,])
    }
    # [nfcstx(m+1)]*[(m+1)x1] = [nfcstx1]
    yhat = newx %*% coef
    # Calculates errors
    loss = abs(sum(ifelse(y[(orig + 1):nT] > yhat, tau, (-1)*(1-tau))))
    # Modifying this part to only print this if specified
    if (print_mdl == 1){
      cat("Losses of out-of-sample forecasts: ", loss, "\n")
    }
  }
  SWfore <- list(coef = coef, yhat = yhat, loss = loss, loadings = M1, 
      DFindex = Dindex, name_vector = nvec)
}

# Decide on the optimal number of vectors.

# (y, x, orig, m, tau) 

#' Function that calculates loss over a given period of time for the diffusion index model
#'
#' @param y - response variables
#' @param x - explanatory variable
#' @param orig - forecast origin
#' @param end - forecasting ending. Note: as the function is currently written on 2/24, this option doesn't do anything.
#' @param m - number of diffusion indices to use
#' @param tau - VaR level
#' @param mod_di - use the modified DI?
#'
#' @return - returns a list of the loss sum and the loss vector
#' @export
#'
#' @examples - loss_calc(pc_df[,1], pc_df[,-1], 757, 1027, 1, 0.01)
loss_calc = function(y, x, orig, m, tau, mod_di = 0, ar_tf = 1, p = 1, print_mdl = 0, model = 1, end = NULL){
  # Extract y_hat values
  if (mod_di == 0){
    di = mod_di(y=y,x=x,orig=orig,m=m, tau=tau, end = end, print_mdl = print_mdl)
  }
  else {
    di = mod_di_wl(y=y,x=x,orig=orig,m=m, tau=tau, ar_tf = ar_tf, p = p, print_mdl = print_mdl, model = model, end = end)
  }
  # mod_di_wl = function (y, x, orig, m, tau, ar_tf = 1, p = 1, print_mdl = 0, model = 1)
  yhat = di$yhat[1:(end-orig)]
  # Calculate the loss
  # Initialize loss vector
  lvec = rep(0,(end-orig))
  # Take the difference
  for (i in 1:(end-orig)){
    # Calculate an indicator variable
    ind = ifelse(y[orig+i] < yhat[i], 1,0)
    # Use indicator in function below
    lvec[i] = (tau - ind)*(y[orig+i] - yhat[i])
  }
  # Add up the losses - change to look at sum of losses. Won't change decision criterion
  sumloss = sum(lvec)
  # sumloss = sum(lvec)/length(lvec)
  return(list(sumloss,lvec))
}


#' Function that selects the optimal number of predictors
#'
#' @param y - response vector
#' @param x - predictor variables
#' @param orig - forecast origin
#' @param end - ending of validation set
#' @param tau - VaR in question
#' @param low_m - low value of m to consider
#' @param high_m - high value of m to consider
#'
#' @return - returns the optimal value of m
#' @export
#'
#' @examples - opt_m(pc_df[,1], pc_df[,-1], 757, 1027, 0.01, low_m =1, high_m  = 5)
opt_m = function(y, x, orig, end = NULL, tau, low_m = 1, high_m, mod_di = 0, ar_tf = 1, p = 1, print_mdl = 0, model = 1, rowname = NULL){
  # Initialize a loss vector
  loss_vec = rep(0,high_m-low_m + 1)
  # Initialize an m vector
  m_vec = seq(low_m, high_m, by = 1)
  # Loop through and populate the loss vector
  for (i in 1:length(loss_vec)){
    loss_vec[i] = quiet(loss_calc(y=y,x=x,orig=orig,end=end, m = m_vec[i], tau = tau, mod_di = mod_di, ar_tf = ar_tf, p = p, print_mdl = print_mdl, model = model))[[1]]
  }
  # Find the minimizer
  opt_m = which.min(loss_vec)
  opt_p = NA
  # Combine into a data frame
  df = as.data.frame(cbind(opt_m, opt_p))
  names(df) <- c("Optimal m", "Optimal p")
  # Assign a rowname
  if (is.null(rowname) == TRUE){
    # Write the row names
    rownames(df) <- c("MV CAViaR")
  }
  else {
    rownames(df) <- rowname
  }
  # Return the loss_vector and the minimzer
  return(list(opt_m, loss_vec, df))
}
#' Function that selects the optimal number of lags
#'
#' @param y - response vector
#' @param x - predictor variables
#' @param orig - forecast origin
#' @param end - ending of validation set
#' @param tau - VaR in question
#' @param low_m - low value of m to consider
#' @param high_m - high value of m to consider
#'
#' @return - returns the optimal value of m
#' @export
#'
#' @examples - opt_mp(y = pc_df[,1], x = pc_df[,-1], orig = 757, end = 1007, tau = 0.01, low_m =1, high_m  = 5, low_p = 1, high_p = 10, ar_tf = 2, mod_di = 1)
opt_mp = function(y, x, orig, end = NULL, tau, low_m = 1, high_m, low_p = 1, high_p, mod_di = 0, ar_tf = 1, print_mdl = 0, model = 1, print_mp = 0, rowname = NULL){
  # Initialize a loss matrix
  loss_mat = matrix(0, high_p-low_p + 1,high_m-low_m + 1)
  # Initialize a p vector
  p_vec = seq(low_p, high_p, by = 1)
  # Loop through and populate the loss vector
  for (i in 1:nrow(loss_mat)){
    loss_mat[i,] = opt_m(y = y, x = x, orig = orig, end = end, tau = tau, low_m = low_m, high_m  = high_m, p = i, mod_di = mod_di, ar_tf = ar_tf, print_mdl = print_mdl, model = model)[[2]]
  }
  # Find the minimizer
  opt_p = which(loss_mat == min(loss_mat), arr.ind = TRUE)[1,1]
  opt_m = which(loss_mat == min(loss_mat), arr.ind = TRUE)[1,2]
  # Print the optimal p and optimal m
  df = as.data.frame(cbind(opt_m, opt_p))
  names(df) <- c("Optimal m", "Optimal p")
  # Assign a rowname
  if (is.null(rowname) == TRUE){
    if (ar_tf == 1){
      # Write the row names
      rownames(df) <- c("MV CAViaR + AR")
    } else if (ar_tf == 2){
      # Write the row names
      rownames(df) <- c("MV CAViaR + SAV")
    } else if (ar_tf == 3){
      # Write the row names
      rownames(df) <- c("MV CAViaR + AS")
    } else {
      rownames(df) <- c("Unknown Model")
    }
  }
  else {
    rownames(df) <- rowname
  }
  # Print the df if the option is turned on
  if (print_mp == 1){
    print(df)
  }
  # Return the loss_vector and the minimzer
  return(list(opt_m, opt_p, loss_mat, df))
}
#' A function that combines optimal values of m and p into a final table
#'
#' @param m1 - the data frame from the "MV CAViaR" run
#' @param m2 - the data frame from the "MV CAViaR + AR" run
#' @param m3 - the data frame from the "MV CAViaR + SAV" run
#' @param m4 - the data frame from the "MV CAViaR + AS" run
#'
#' @return - a nicely formatted table
#' @export
#'
#' @examples - pretty_pm(opt_pred_nl[[3]], opt_pm_m1[[4]], opt_pm_m2[[4]], opt_pm_m3[[4]])
pretty_pm = function(m1, m2, m3, m4){
  # Merge the individual data frames
  pm_pretty_df = rbind(m1, m2, m3, m4)
  # Format nicely
  pm_pretty_df %>% kable(caption = "Optimal Number of Diffusion Indices (m) and Lags (p) for Different Models", digits = 0) %>% kable_styling("striped", full_width = F) %>% kable_styling() %>% footnote(general = "The MV CAViaR model doesn't have an optimal value for p because there are no AR lags in the model"
  )
  
}

#' Here is a function that runs the univariate CAViaR model 4 times
#'
#' @param df - the percent change data frame to consider
#' @param nfcst - number of forecasts to run
#' @param tau - the VaR level to consider
#' @param no_run - specifies if any models should not be run
#'
#' @return - a list of the 4 univariate model forecasts
#' @export
#'
#' @examples - aceg = gen_uv_test(pc_df, 1, 0.05, no_run = c(1,1,0,1))
gen_uv_test = function(df, nfcst, tau, no_run = c(0,0,0,0)){
  # model type (1 - SAV, 2 - AS, 3 - GARCH, 4 - ADAPTIVE) 
  # Initialize a list  
  out_list = list()
  # Run the four models - model 1; SAV
  if (no_run[1] == 0){
    uvcav_1 = rolling_predictions(df[,1], range_data = (1:length(df[,1])), nfcst = nfcst, model = 1, G = 10, col = 1, level = tau)
  }
  # Add a filler if there's no entry
  else {
    uvcav_1 = 0
  }
  # Model 2 - AS
  if (no_run[2] == 0){
    uvcav_2 = rolling_predictions(df[,1], range_data = (1:length(df[,1])), nfcst = nfcst, model = 2, G = 10, col = 1, level = tau)
  }
  else {
    uvcav_2 = 0
  }
  # Model 3 - GARCH
  if (no_run[3] == 0){
    uvcav_3 = rolling_predictions(df[,1], range_data = (1:length(df[,1])), nfcst = nfcst, model = 3, G = 10, col = 1, level = tau)
  }
  else {
    uvcav_3 = 0
  }
  # Model 4 - Adaptive
  if (no_run[4] == 0){
    uvcav_4 = rolling_predictions(df[,1], range_data = (1:length(df[,1])), nfcst = nfcst, model = 4, G = 10, col = 1, level = tau)
  }
  else {
    uvcav_4 = 0
  }
  # Export the data as a list
  return(list(uvcav_1, uvcav_2, uvcav_3, uvcav_4))
}

#' Function to plot the data which we generate in previous functions
#'
#' @param plot_matrix - matrix with the data to plot
#' @param norm_value - what to subtact from the data to make it on a percentage change basis. Default is 100.
#'
#' @return 
#' @export - a plot of the data by diffusion index number
#'
#' @examples = plt_data(plot_mtx[[1]]), abc = plt_data(plot_mat, tau = 0.01)
plt_data = function(plot_matrix, tau, resp_var, ntest){
  # Establish a maximum and minimum value
  max_val = max(plot_matrix[,1:ncol(plot_matrix)])
  min_val = min(plot_matrix[,1:ncol(plot_matrix)])
  # Calculate inital and ending time value
  start = index(plot_matrix)[1]
  end = index(plot_matrix)[nrow(plot_matrix)]
  ind_vals = index(plot_matrix) - start
  # Create an initial plot and add lines
    for (i in 1:ncol(plot_matrix)){
      if (i == 1){
        # 4/2/2020 - fixing the index
        plot.ts(ind_vals,plot_matrix[,i], type = "l", xlab = paste("Days Since", as.Date(start)), ylab = "Percent Change in PG", ylim = c(min_val,max_val), lwd = 1, main = paste("Predicting", resp_var, "Returns from", as.Date(start), "to", as.Date(end)), sub = paste("The VaR Level is ", 100*tau, "%", "; There are ", ntest, " Trading Days Plotted Above", sep = ""))
        # plot.ts(index(plot_matrix), plot_matrix[,i], type = "l", xlab = "Trading Days", ylab = "Percent Change in PG", ylim = c(min_val,max_val), lwd = 1, main = "Predicting PG Returns Over Last 250 Trading Days in 2008", sub = paste("The VaR Level is ", 100*tau, "%", sep = ""))
    } else if(i %in% seq(2,8,1)) {
        lines(ind_vals,plot_matrix[,i], col = i-1, lty = 2)
    } else {
        lines(ind_vals,plot_matrix[,i], col = i-1, lty = 2, lwd = 2)
    }
    }
  # Define a sequence for plotting
  plot_seq = seq(1, ncol(plot_matrix))
  legend("topleft", legend = c(colnames(plot_matrix)), col = plot_seq, lty = c(1, rep(2, 7), rep(3, ifelse(ncol(plot_matrix)-8 <= 0, 0, ncol(plot_matrix)-8))), lwd = c(1, rep(1, 7), rep(2, ifelse(ncol(plot_matrix)-8 <= 0, 0, ncol(plot_matrix)-8))))
  # Add a line for 0
  # abline(h = 0, col = "black", lty = 2)
}


#' A function to calculate losses based on the test sample
#'
#' @param true_vec - the true vector of returns
#' @param pred_vec - the predicted vector from the model runs
#' @param tau - VaR level. Must match what the model used
#'
#' @return - total losses and the entire loss vector
#' @export
#'
#' @examples
loss_test = function(true_vec, pred_vec, tau){
  # Initialize a loss vector
  lvec = rep(0, length(true_vec))
  # Initialize a break vector to see when VaR is broken
  bvec = rep(0, length(true_vec))
  for (i in 1:length(true_vec)){
    # Calculate an indicator variable
    bvec[i] = ifelse(true_vec[i] < pred_vec[i], 1,0)
    # Use indicator in function below
    lvec[i] = (tau - bvec[i])*(true_vec[i] - pred_vec[i])
  }
  # Add up the losses
  # sumloss = sum(lvec)/length(lvec)
  sumloss = sum(lvec)
  # Add up the VaR breakage
  varbreak = sum(bvec)/length(bvec)
  return(list(sumloss,lvec, varbreak, bvec))
}
#' A function to calculate losses based on the plot matrix
#'
#' @param data_mat - a matrix of forecasted VaR values, with the true value in the first column
#' @param tau - VaR level. Must match what the model used
#'
#' @return - a list of four items. 
#' 1 = a vector of the losses of all models. 
#' 2 = a vector showing the percentage of VaR breaks by model
#' 3 = the loss matrix
#' 4 = the break matrix
#' @export
#'
#' @examples
gen_loss_test = function(data_mat, tau){
  # Initialize loss and break matrices
  lmat = bmat = matrix(0, nrow = nrow(data_mat), ncol = ncol(data_mat)-1)
  # bvec = rep(0, length(true_vec))
  # Populate the matrices
  for (i in 1:nrow(lmat)){
    for (j in 1:(ncol(lmat))){
      # Calculate an indicator variable
      bmat[i,j] = ifelse(data_mat[i,1] < data_mat[i,j+1], 1,0)
      # Use indicator in function below
      lmat[i,j] = (tau - bmat[i,j])*(data_mat[i,1] - data_mat[i,j+1])
    }    
  }
  # Add up the losses
  sumloss = colSums(lmat)
  # Add up the VaR breakage
  varbreak = colSums(bmat)/nrow(bmat)
  return(list(sumloss, varbreak, lmat, bmat))
}
#' A function to make a nice comparison of losses
#'
#' @param data_mat - input data matrix used in the calculation of losses
#' @param loss_list - a list of the losses calculated from the CAViaR function
#' @param tau - the risk level used
#' @param ntest - the number of test points
#'
#' @return
#' @export - returns a nicely formatted table
#'
#' @examples - pretty_tables(plot_mat, l_list, tau = 0.01)
pretty_tables = function(data_mat, loss_list, tau, ntest){
  # Combine into a data frame
  df = as.data.frame(rbind(loss_list[[1]], loss_list[[2]]))
  # Calculate inital and ending time value
  start = index(data_mat)[1]
  end = index(data_mat)[nrow(data_mat)]
  # Add row/column names
  colnames(df) <- colnames(data_mat[,-1])
  rownames(df) <- c("Losses", "VaR Breaks (%)")
  # Convert to a table
  df %>% kable(caption = paste("Comparison of VaR Methods for a ", tau*100, "% VaR", sep = ""), digits = 3) %>% kable_styling("striped", full_width = F) %>% kable_styling() %>% footnote(general = paste("Calculated using", ntest, "trading days from", as.Date(start), "to", as.Date(end)))
}
#' A dressed up version of the export function
#'
#' @param var_file - file to export
#' @param path - filepath
#' @param filename - name of the file, ending with .CSV
#'
#' @return
#' @export - exported CSV file
#'
#' @examples - exp_func(var_file = var_1pc_2016_usetf[[1]], path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "TEST.csv")
exp_func = function(var_file, path, filename){
  # Write a zoo
  write.zoo(var_file, paste0(path, filename), quote = FALSE, sep = ",")
}

# exp_func(var_file = var_1pc_2016_usetf[[1]], path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "TEST.csv")

#' This is the "master" function where we'll evaluate the importance of the VaR model over several time periods
#'
#' @param symbol_list - a list of symbols to feed into the model 
#' @param resp_var - the response variable
#' @param compl_case - should the model require complete cases? Default value is 1.
#' @param adj_close - use adjusted close price for the predictors? Default value is 1.
#' @param resp_adj_close - use adjusted close price for the response? Default value is 1.
#' @param start_date - start date to pull data from
#' @param end_date - end date to pull data from
#' @param nval - number of validation points to use
#' @param ntest - number of test points to use
#' @param tau - VaR level to use
#' @param low_m - low number of predictors to test
#' @param high_m  - low number of predictors to test
#' @param uv_list - a list of a pre-run univariate model. If a data frame is not provided, the lengthy uv model will run 
#' @param no_run - things not to run in the model
#' @param low_p - low value for number of lags
#' @param high_p - high value for number of lags
#' @param na_interp - should the function interpolate NA's
#' @param print_mdl - print the model summaries?
#' @param print_mp - print the optimal values for p and m
#' @param lag_pred - do you want to lag the m predictors (default is 1; strongly recommended)
#' @param rowname - what to name the rows of the nice p and m matrix
#' @param export_csv - do you want to export a CSV? Default is 1.
#' @param path - path to export the CSV
#' @param filename - what to name the CSV
#'
#' @return - a list of the plot matrix, a plot, a list with losses, and a table
#' @export - a plot and a table
#'
#' @examples - cav_simul(c("DIS", "GE", "IBM", "MMM", "XOM"), resp_var = "PG", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 5, tau = 0.01, uv_list = uv_cav_list)
cav_simul = function(symbol_list, resp_var, compl_case = 1, adj_close = 1, resp_adj_close = 1, start_date = "1900-01-01", end_date = Sys.Date(), nval, ntest, tau, low_m = 1, high_m, low_p = 1, high_p, uv_list = NULL, no_run = c(0,0,0,0), na_interp = TRUE, print_mdl = 0, print_mp = 0, lag_pred = 1, rowname = NULL, export_csv = 1, path, filename){
  # Select data parameters, pull the data, and percent change the data
  df = diff_index_df(symbol_list = symbol_list, resp_var = resp_var, compl_case = compl_case, adj_close = adj_close, resp_adj_close = resp_adj_close, start_date = start_date, end_date = end_date, lag_pred = lag_pred)
  # Take the percent change of the data
  pc_df = pc_diff_index(df)
  # Extract the legnth of the data frame
  nr = test_end = nrow(pc_df)
  # Calculate the start of the val period, the end of the val period, and the beginning and end of test period
  test_orig = test_end - ntest
  val_end = test_orig
  val_orig = test_orig - nval
  # Test for the optimal number of parameters
  opt_pred_nl = opt_m(y = pc_df[,1], x = pc_df[,-1], orig = val_orig, end = val_end, tau = tau, low_m = low_m, high_m = high_m, rowname = rowname)
  opt_pred_art1 = opt_mp(y = pc_df[,1], x = pc_df[,-1], orig = val_orig, end = val_end, tau = tau, low_m = low_m, high_m = high_m, low_p = low_p, high_p = high_p,mod_di = 1, ar_tf = 1, print_mdl = print_mdl, print_mp = print_mp, rowname = rowname)
  opt_pred_art2 = opt_mp(y = pc_df[,1], x = pc_df[,-1], orig = val_orig, end = val_end, tau = tau, low_m = low_m, high_m = high_m, low_p = low_p, high_p = high_p,mod_di = 1, ar_tf = 2, print_mdl = print_mdl, print_mp = print_mp, rowname = rowname)
  opt_pred_art3 = opt_mp(y = pc_df[,1], x = pc_df[,-1], orig = val_orig, end = val_end, tau = tau, low_m = low_m, high_m = high_m, low_p = low_p, high_p = high_p, mod_di = 1, ar_tf = 3, print_mdl = print_mdl, print_mp = print_mp, rowname = rowname)
  # gen_uv_test(pc_df, 1, 0.05, no_run = c(1,1,0,1))
  # Use the above forecasts to input into the above
  mv_fcst = mod_di(pc_df[,1], pc_df[,-1], orig = test_orig, m = opt_pred_nl[[1]], tau = tau, print_mdl = print_mdl)
  mv_fcst_art1 = mod_di_wl(pc_df[,1], pc_df[,-1], orig = test_orig, m = opt_pred_art1[[1]], p = opt_pred_art1[[2]], tau = tau, ar_tf = 1, print_mdl = print_mdl)
  mv_fcst_art2 = mod_di_wl(pc_df[,1], pc_df[,-1], orig = test_orig, m = opt_pred_art2[[1]], p = opt_pred_art2[[2]], tau = tau, ar_tf = 2, print_mdl = print_mdl)
  mv_fcst_art3 = mod_di_wl(pc_df[,1], pc_df[,-1], orig = test_orig, m = opt_pred_art3[[1]], p = opt_pred_art3[[2]], tau = tau, ar_tf = 3, print_mdl = print_mdl)
  # Calculate the number of predictions
  if (is.null(uv_list) == TRUE){
    # Print a warning
    print("WARNING: Not supplying an input data frame will require this function to run for a significant amount of time (1hr+)")
    # Call the function
    # gen_uv_test = function(df, nfcst, tau, no_run = c(0,0,0,0)){
    # print(head(pc_df))
    uv_list = gen_uv_test(df = pc_df, nfcst = ntest, tau = tau, no_run = no_run)
    # Add to a data frame
    # Incorporate the rolling predictions function results here
  plot_mat = cbind(pc_df[(test_orig+1):nrow(pc_df),1], mv_fcst$yhat[1:ntest], mv_fcst_art1$yhat[1:ntest], mv_fcst_art2$yhat[1:ntest], mv_fcst_art3$yhat[1:ntest], uv_list[[1]][(test_orig+1):test_end]*(-1), uv_list[[2]][(test_orig+1):test_end]*(-1), uv_list[[3]][(test_orig+1):test_end]*(-1), uv_list[[4]][(test_orig+1):test_end]*(-1))
  } else {
    # Assign the columns of the data frame
    # head(var_5pc_2010_usetf[[1]][,6:9])
    # model type (1 - SAV, 2 - AS, 3 - GARCH, 4 - ADAPTIVE) 
    # test_df = head(var_5pc_2010_usetf[[1]][,6:9])
    # test_df$SAV
    # test_df$`Abs. Slope`
    # test_df$`Ind. GARCH`
    # test_df$Adaptive
    plot_mat = cbind(pc_df[(test_orig+1):nrow(pc_df),1], mv_fcst$yhat[1:ntest], mv_fcst_art1$yhat[1:ntest], mv_fcst_art2$yhat[1:ntest], mv_fcst_art3$yhat[1:ntest], uv_list$SAV, uv_list$`Abs. Slope`, uv_list$`Ind. GARCH`, uv_list$Adaptive)
  }
  # Count the NAs and print a warning
  print(paste("NOTE: There are ", sum(is.na(plot_mat)), " NA(s) in the dataset", sep = ""))
  # Linearly interpolate the NAs
  if (na_interp == TRUE){
    # Assign the plot matrix to a new value
    plot_mat_na <- plot_mat
    # Print a warning
    print("WARNING: There were missing values in the plot matrix.")
    # Interpolate the NA's
    for (i in 1:ncol(plot_mat_na)){
      # Interpolate the data
      plot_mat[,i] <- na.approx(plot_mat_na[,i])
    }
  }
  # model type (1 - SAV, 2 - AS, 3 - GARCH, 4 - ADAPTIVE)
  # Add descriptive titles onto the plot_mat
  colnames(plot_mat) <- c(resp_var, "MV CAViaR", "MV CAViaR + AR", "MV CAViaR + SAV", "MV CAViaR + AS", "SAV", "Abs. Slope", "Ind. GARCH", "Adaptive")
  # colnames(plot_mat) <- c("SPY", "MV CAViaR", "MV CAViaR + AR", "MV CAViaR + SAV", "MV CAViaR + AS", "SAV", "Abs. Slope", "Ind. GARCH", "Adaptive")
  # Plot everything
  plot = plt_data(plot_mat, tau = tau, resp_var = resp_var, ntest = ntest)
  # Calculate losses
  l_list = gen_loss_test(plot_mat, tau = tau)
  # Put into tables
  tables = pretty_tables(plot_mat, l_list, tau = tau, ntest = ntest)
  # Run the function for optimal p and m
  pm_table = pretty_pm(opt_pred_nl[[3]], opt_pred_art1[[4]], opt_pred_art2[[4]], opt_pred_art3[[4]])
  # Export the matrix
  if (export_csv == 1){
    exp_func(var_file = plot_mat, path, filename)
  }
  # Print the tables and the plot
  print(plot)
  print(tables)
  print(pm_table)
  return(list(plot_mat, plot, l_list, tables, plot_mat_na, pm_table))
}
#' A function to input the VaR files, plot them and generate tables
#'
#' @param file_path - file path to use
#' @param filename - name of the file
#' @param tau - quantile to use
#' @param resp_var - response variable to use in the plot
#' @param ntest - number of test points
#' @param cn_input - column name inputs
#'
#' @return - a list of the xts file, the plot, the loss list, and tables
#' @export - a plot and tables
#'
#' @examples - test = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_1pc_2008_us_etf.csv", 0.01)
var_input_disp = function(file_path, filename, tau, resp_var = "SPY", ntest = 250, cn_input = c("SPY", "M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8")){
# var_input_disp = function(file_path, filename, tau, resp_var = "SPY", ntest = 250, cn_input = c("SPY", "MV CAViaR", "MV CAViaR + AR", "MV CAViaR + SAV", "MV CAViaR + AS", "SAV", "Abs. Slope", "Ind. GARCH", "Adaptive")){
  # Import data
  plot_mat = read.csv(paste0(file_path,filename), sep = ",", header = T, stringsAsFactors = FALSE)
  # Fix date format
  plot_mat$Index = as.Date(plot_mat$Index)
  # Convert to an xts
  plot_mat = xts(plot_mat[,-1], order.by = plot_mat[,1])
  # Fix column names
  colnames(plot_mat) <- cn_input
  # Plot everything
  plt_data(plot_mat, tau = tau, resp_var = resp_var, ntest = ntest)
  # plot = plt_data(plot_mat, tau = tau, resp_var = resp_var, ntest = ntest)
  # Calculate losses
  l_list = gen_loss_test(plot_mat, tau = tau)
  # Put into tables
  pm_table = pretty_tables(plot_mat, l_list, tau = tau, ntest = ntest)
  # Print the tables and the plot
  # print(plot)
  # print(tables)
  print(pm_table)
  # Return the xts, the plot, the loss list, and the tables
  return(list(plot_mat))
  # return(list(plot_mat, plot, l_list, tables))
}
# return(list(plot_mat, plot, l_list, tables))

# Call the above function
v1_2008_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_1pc_2008_us_etf.csv", 0.01)
# pretty_tables(v1_2008_usetf[[1]], v1_2008_usetf[[3]], tau = 0.01, ntest = 250)
v5_2008_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_5pc_2008_us_etf.csv", 0.05)
v10_2008_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_10pc_2008_us_etf.csv", 0.10)

# head(v10_2008_usetf[[1]])
# v10_2008_usetf
# 1%, 5%, 10% VaR - 2008 - 1st set of predictors
# var_1pc_2008_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_1pc_2008_us_etf.csv", uv_list = var_1pc_2008_usetf[[1]])

# var_5pc_2008_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_5pc_2008_us_etf.csv", uv_list = var_5pc_2008_usetf[[1]])

# var_10pc_2008_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_10pc_2008_us_etf.csv", uv_list = var_10pc_2008_usetf[[1]])
# Call the above function
v1_2008_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_1pc_2008_glob_etf.csv", 0.01)
v5_2008_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_5pc_2008_glob_etf.csv", 0.05)
v10_2008_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_10pc_2008_glob_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2008 - 2nd set of predictors
# var_1pc_2008_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_1pc_2008_glob_etf.csv", uv_list = var_1pc_2008_usetf[[1]])

# var_5pc_2008_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_5pc_2008_glob_etf.csv", uv_list = var_5pc_2008_usetf[[1]])

# var_10pc_2008_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_10pc_2008_glob_etf.csv", uv_list = var_10pc_2008_usetf[[1]])
# 1%, 5%, 10% VaR - 2008 - 3rd set of predictors
# var_1pc_2008_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_1pc_2008_comm_etf.csv", uv_list = var_1pc_2008_usetf[[1]])

# var_5pc_2008_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_5pc_2008_comm_etf.csv", uv_list = var_5pc_2008_usetf[[1]])

# var_10pc_2008_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_10pc_2008_comm_etf.csv", uv_list = var_10pc_2008_usetf[[1]])
# Call the above function
v1_2008_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_1pc_2008_bond_etf.csv", 0.01)
v5_2008_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_5pc_2008_bond_etf.csv", 0.05)
v10_2008_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_10pc_2008_bond_etf.csv", 0.10)
# iShares 1-3 Year Treasury Bond Fund (SHY)
# iShares 7-10 Year Treasury Bond Fund (IEF)
# iShares 20+ Year Treasury Bond Fund (TLT)
# iShares iBoxx $ Investment Grade Corporate Bond ETF (LQD)

# 1%, 5%, 10% VaR - 2008 - 4th set of predictors
# var_1pc_2008_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_1pc_2008_bond_etf.csv", uv_list = var_1pc_2008_usetf[[1]])

# var_5pc_2008_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_5pc_2008_bond_etf.csv", uv_list = var_5pc_2008_usetf[[1]])

# var_10pc_2008_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_10pc_2008_bond_etf.csv", uv_list = var_10pc_2008_usetf[[1]])
# Call the above function
v1_2008_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_1pc_2008_all_etf.csv", 0.01)
v5_2008_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_5pc_2008_all_etf.csv", 0.05)
v10_2008_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_10pc_2008_all_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2008 - 5th set of predictors
# var_1pc_2008_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_1pc_2008_all_etf.csv", uv_list = var_1pc_2008_usetf[[1]])

# var_5pc_2008_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_5pc_2008_all_etf.csv", uv_list = var_5pc_2008_usetf[[1]])

# var_10pc_2008_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2004-01-01", end_date = "2008-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_10pc_2008_all_etf.csv", uv_list = var_10pc_2008_usetf[[1]])
# Call the above function
v1_2010_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_1pc_2010_us_etf.csv", 0.01)
v5_2010_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_5pc_2010_us_etf.csv", 0.05)
v10_2010_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_10pc_2010_us_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2010 - 1st set of predictors
# var_1pc_2010_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_1pc_2010_us_etf.csv", uv_list = var_1pc_2010_usetf[[1]])

# var_5pc_2010_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_5pc_2010_us_etf.csv", uv_list = var_5pc_2010_usetf[[1]])

# var_10pc_2010_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_10pc_2010_us_etf.csv", uv_list = var_10pc_2010_usetf[[1]])
# Call the above function
v1_2010_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_1pc_2010_glob_etf.csv", 0.01)
v5_2010_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_5pc_2010_glob_etf.csv", 0.05)
v10_2010_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_10pc_2010_glob_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2010 - 2nd set of predictors
# var_1pc_2010_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_1pc_2010_glob_etf.csv", uv_list = var_1pc_2010_usetf[[1]])

# var_5pc_2010_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_5pc_2010_glob_etf.csv", uv_list = var_5pc_2010_usetf[[1]])

# var_10pc_2010_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_10pc_2010_glob_etf.csv", uv_list = var_10pc_2010_usetf[[1]])
# 1%, 5%, 10% VaR - 2010 - 3rd set of predictors
# var_1pc_2010_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_1pc_2010_comm_etf.csv", uv_list = var_1pc_2010_usetf[[1]])

# var_5pc_2010_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_5pc_2010_comm_etf.csv", uv_list = var_5pc_2010_usetf[[1]])

# var_10pc_2010_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_10pc_2010_comm_etf.csv", uv_list = var_10pc_2010_usetf[[1]])
# Call the above function
v1_2010_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_1pc_2010_bond_etf.csv", 0.01)
v5_2010_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_5pc_2010_bond_etf.csv", 0.05)
v10_2010_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_10pc_2010_bond_etf.csv", 0.10)
# iShares 1-3 Year Treasury Bond Fund (SHY)
# iShares 7-10 Year Treasury Bond Fund (IEF)
# iShares 20+ Year Treasury Bond Fund (TLT)
# iShares iBoxx $ Investment Grade Corporate Bond ETF (LQD)

# 1%, 5%, 10% VaR - 2010 - 4th set of predictors
# var_1pc_2010_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_1pc_2010_bond_etf.csv", uv_list = var_1pc_2010_usetf[[1]])

# var_5pc_2010_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_5pc_2010_bond_etf.csv", uv_list = var_5pc_2010_usetf[[1]])

# var_10pc_2010_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_10pc_2010_bond_etf.csv", uv_list = var_10pc_2010_usetf[[1]])
# Call the above function
v1_2010_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_1pc_2010_all_etf.csv", 0.01)
v5_2010_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_5pc_2010_all_etf.csv", 0.05)
v10_2010_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_10pc_2010_all_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2010 - 5th set of predictors
# var_1pc_2010_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_1pc_2010_all_etf.csv", uv_list = var_1pc_2010_usetf[[1]])

# var_5pc_2010_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_5pc_2010_all_etf.csv", uv_list = var_5pc_2010_usetf[[1]])

# var_10pc_2010_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2006-01-01", end_date = "2010-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_10pc_2010_all_etf.csv", uv_list = var_10pc_2010_usetf[[1]])
# Call the above function
v1_2014_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_1pc_2014_us_etf.csv", 0.01)
v5_2014_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_5pc_2014_us_etf.csv", 0.05)
v10_2014_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_10pc_2014_us_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2014 - 1st set of predictors
# var_1pc_2014_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_1pc_2014_us_etf.csv", uv_list = var_1pc_2014_usetf[[1]])

# var_5pc_2014_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_5pc_2014_us_etf.csv", uv_list = var_5pc_2014_usetf[[1]])

# var_10pc_2014_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_10pc_2014_us_etf.csv", uv_list = var_10pc_2014_usetf[[1]])
# Call the above function
v1_2014_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_1pc_2014_glob_etf.csv", 0.01)
v5_2014_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_5pc_2014_glob_etf.csv", 0.05)
v10_2014_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_10pc_2014_glob_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2014 - 2nd set of predictors
# var_1pc_2014_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_1pc_2014_glob_etf.csv", uv_list = var_1pc_2014_usetf[[1]])

# var_5pc_2014_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_5pc_2014_glob_etf.csv", uv_list = var_5pc_2014_usetf[[1]])

# var_10pc_2014_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_10pc_2014_glob_etf.csv", uv_list = var_10pc_2014_usetf[[1]])
# 1%, 5%, 10% VaR - 2014 - 3rd set of predictors
# var_1pc_2014_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_1pc_2014_comm_etf.csv", uv_list = var_1pc_2014_usetf[[1]])

# var_5pc_2014_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_5pc_2014_comm_etf.csv", uv_list = var_5pc_2014_usetf[[1]])

# var_10pc_2014_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_10pc_2014_comm_etf.csv", uv_list = var_10pc_2014_usetf[[1]])
# Call the above function
v1_2014_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_1pc_2014_bond_etf.csv", 0.01)
v5_2014_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_5pc_2014_bond_etf.csv", 0.05)
v10_2014_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_10pc_2014_bond_etf.csv", 0.10)
# iShares 1-3 Year Treasury Bond Fund (SHY)
# iShares 7-10 Year Treasury Bond Fund (IEF)
# iShares 20+ Year Treasury Bond Fund (TLT)
# iShares iBoxx $ Investment Grade Corporate Bond ETF (LQD)

# 1%, 5%, 10% VaR - 2014 - 4th set of predictors
# var_1pc_2014_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_1pc_2014_bond_etf.csv", uv_list = var_1pc_2014_usetf[[1]])

# var_5pc_2014_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_5pc_2014_bond_etf.csv", uv_list = var_5pc_2014_usetf[[1]])

# var_10pc_2014_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_10pc_2014_bond_etf.csv", uv_list = var_10pc_2014_usetf[[1]])
# Call the above function
v1_2014_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_1pc_2014_all_etf.csv", 0.01)
v5_2014_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_5pc_2014_all_etf.csv", 0.05)
v10_2014_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_10pc_2014_all_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2014 - 5th set of predictors
# var_1pc_2014_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_1pc_2014_all_etf.csv", uv_list = var_1pc_2014_usetf[[1]])

# var_5pc_2014_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_5pc_2014_all_etf.csv", uv_list = var_5pc_2014_usetf[[1]])

# var_10pc_2014_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2010-01-01", end_date = "2014-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_10pc_2014_all_etf.csv", uv_list = var_10pc_2014_usetf[[1]])
# Call the above function
v1_2016_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_1pc_2016_us_etf.csv", 0.01)
v5_2016_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_5pc_2016_us_etf.csv", 0.05)
v10_2016_usetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/","var_10pc_2016_us_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2016 - 1st set of predictors
# var_1pc_2016_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_1pc_2016_us_etf.csv", uv_list = var_1pc_2016_usetf[[1]])

# var_5pc_2016_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_5pc_2016_us_etf.csv", uv_list = var_5pc_2016_usetf[[1]])

# var_10pc_2016_us_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 9, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_US_ETF_runs/", filename = "var_10pc_2016_us_etf.csv", uv_list = var_10pc_2016_usetf[[1]])
# Call the above function
v1_2016_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_1pc_2016_glob_etf.csv", 0.01)
v5_2016_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_5pc_2016_glob_etf.csv", 0.05)
v10_2016_globetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/","var_10pc_2016_glob_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2016 - 2nd set of predictors
# var_1pc_2016_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_1pc_2016_glob_etf.csv", uv_list = var_1pc_2016_usetf[[1]])

# var_5pc_2016_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_5pc_2016_glob_etf.csv", uv_list = var_5pc_2016_usetf[[1]])

# var_10pc_2016_glob_etf = cav_simul(c("JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 10, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_glob_ETF_runs/", filename = "var_10pc_2016_glob_etf.csv", uv_list = var_10pc_2016_usetf[[1]])
# 1%, 5%, 10% VaR - 2016 - 3rd set of predictors
# var_1pc_2016_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_1pc_2016_comm_etf.csv", uv_list = var_1pc_2016_usetf[[1]])

# var_5pc_2016_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_5pc_2016_comm_etf.csv", uv_list = var_5pc_2016_usetf[[1]])

# var_10pc_2016_comm_etf = cav_simul(c("DBA", "DBC", "DBE", "DBB"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_comm_ETF_runs/", filename = "var_10pc_2016_comm_etf.csv", uv_list = var_10pc_2016_usetf[[1]])
# Call the above function
v1_2016_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_1pc_2016_bond_etf.csv", 0.01)
v5_2016_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_5pc_2016_bond_etf.csv", 0.05)
v10_2016_bondetf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/","var_10pc_2016_bond_etf.csv", 0.10)
# iShares 1-3 Year Treasury Bond Fund (SHY)
# iShares 7-10 Year Treasury Bond Fund (IEF)
# iShares 20+ Year Treasury Bond Fund (TLT)
# iShares iBoxx $ Investment Grade Corporate Bond ETF (LQD)

# 1%, 5%, 10% VaR - 2016 - 4th set of predictors
# var_1pc_2016_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_1pc_2016_bond_etf.csv", uv_list = var_1pc_2016_usetf[[1]])

# var_5pc_2016_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_5pc_2016_bond_etf.csv", uv_list = var_5pc_2016_usetf[[1]])

# var_10pc_2016_bond_etf = cav_simul(c("SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 4, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_bond_ETF_runs/", filename = "var_10pc_2016_bond_etf.csv", uv_list = var_10pc_2016_usetf[[1]])
# Call the above function
v1_2016_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_1pc_2016_all_etf.csv", 0.01)
v5_2016_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_5pc_2016_all_etf.csv", 0.05)
v10_2016_alletf = var_input_disp("/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/","var_10pc_2016_all_etf.csv", 0.10)
# 1%, 5%, 10% VaR - 2016 - 5th set of predictors
# var_1pc_2016_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10, tau = 0.01, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_1pc_2016_all_etf.csv", uv_list = var_1pc_2016_usetf[[1]])

# var_5pc_2016_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.05, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_5pc_2016_all_etf.csv", uv_list = var_5pc_2016_usetf[[1]])

# var_10pc_2016_all_etf = cav_simul(c("XLU", "XLP", "XLV", "XLK", "XLY", "XLI", "XLF", "XLB", "XLE", "JXI", "KXI", "IXJ", "IXP", "IXN", "RXI", "EXI", "IXG", "MXI", "IXC", "SHY", "IEF", "TLT", "LQD"), resp_var = "SPY", start_date = "2012-01-01", end_date = "2016-12-31", nval = 250, ntest = 250, low_m = 1, high_m = 23, low_p = 1, high_p = 10,tau = 0.10, print_mdl = 1, print_mp = 1, path = "/Users/stevenmoen/Documents/GitHub/CAViaR_MS_thesis/Data_Export/SPY_all_ETF_runs/", filename = "var_10pc_2016_all_etf.csv", uv_list = var_10pc_2016_usetf[[1]])

Literature Cited